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HIGH ORDER FINITE DIFFERENCE METHODS. MULTIDIMENSIONAL LINEAR 
PROBLEMS AND CURVILINEAR COORDINATES 


JAN NORDSTROM* AND MARK H. CARPENTER* 

Abstract. Boundary and interface conditions are derived for high order finite difference methods applied 
to multidimensional linear problems in curvilinear coordinates. The boundary and interface conditions lead 
to conservative schemes and strict and strong stability provided that certain metric conditions are met. 

Key words, high-order finite-difference, numerical stability, interface conditions, summation-by-parls. 
variable coefficient 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. Phenomena that require an accurate description of high frequency variation in space 
for long times occur in many important applications such as electromagnetics, acoustics (all cases of wave 
propagation), and direct simulation of turbulent and transitional flow; see for example [l]-[6]. Strictly stable 
high order finite difference methods are well suited for these types of problems (see [7]- [16]) because they 
guarantee accurate results with bounded error growth in time for realistic meshes. 

Most of the development for these types of methods has considered constant coefficient problems on a 
Cartesian mesh. In [17], [18] stable and conservative boundary and interface conditions were derived for 
the one-dimensional constant coefficient Euler and Navier-Stokes equations on multiple domains. A similar 
technique was used in [19], [20], and [21] for Chebyshev spectral methods. 

In this paper we extend the constant coefficient analysis in [17], [18] to scalar multidimensional linear 
problems in curvilinear coordinates including block interfaces. Related previous work includes investigations 
of the metric derivatives in non-smooth meshes (see [22], [23]) and the treatment of parabolic and hyperbolic 
systems in curvilinear coordinates on a single domain [14]. 

The rest of this paper will proceed as follows. Section 2, will give some basic definitions. Section 3 
presents the ID difference operators that form the basis of the multidimensional approximation treatment. 
Section 4 defines the linear model problem and discusses we 11- posed ness. Section 5 provides an investigation 
of the discrete problem. Section 6 illustrates numerical experiments and in Section 7 we summarize and 
draw conclusions. 

2. Definitions. Consider the linear initial boundary value problem 


w t ~ P(jL .t)iv -f 8F(x,i) 

, x G ft 

J > 0 , 

w = Sf(r) 

, r G fl 

./ = (), 

Lrw = 6<l(t) 

, .r G F 

J > 0 , 


where P is the differential operator, Lc is the boundary operator, Q is the domain and I is the domain 
boundary. The forcing function tfF, the initial function Sf, and the boundary data 8g are the data of the 
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problem; w denotes the difference between a solution with data /, F, g and one with data / + Sf, F + SF, 
9 + Sg. 

The semi-discrete version of (2.1) is 


(U’j)t = Q{ 2 j,t)wj +6Fj(t) 

, xj € ft 

J > 0, 

«'j = s fj 

, x.j e ft 

,t = 0, 

L D Wj = Sg(t) 

- * j € r 

J > 0, 


where Q is the difference operator approximating the differential operator P, SFj is the forcing function, Sfj 
the initial function, L& the discrete boundary operator where numerical boundary conditions are included, 
and Sg the boundary data. It is assumed that (2.2) is a consistent approximation of (2.1). 

2.1. Well-posedness and stability. There are many concepts of we 11- posed ness and stability; see 
[24]. Here we consider the following definitions. 

Definition 1. Problem (2.1) is strongly well posed if the solution w is unique, exists , and satisfies 

(2-3) IMIn+ / IHIr* < A'^IIMIn + /Wlln + \M\l)dt}, 

Jo Jo 

where I\ c and i} c may not depend on SF , Sf , Sg, and || * ||^ and || * ||r are suitable continuous norms. 

Definition 2. Problem (2.2) is strongly stable if for a sufficiently fine mesh the solution wj satisfies 

(2.4) IHIn + f IIHIr<« < W* 1 WfUh + All^llh + ll^llf )*}, 

Jo Jo 

where I\ d and g d may not depend on SFj , Sfj, Sg , and || • ||n and || * ||r art suitable discrete norms. 

Definition 3. The approximation (2.2) of (2.1) is strictly stable if the analytical and discrete growth 
rates (see (2.3) and (2.4)) satisfy 

(2.5) 7) d < T] c -f (?(Ax), 


where Ax is the mesh size. 


2.2. Linear algebra relations. For later reference we define some useful matrix operations; see [25]. 
DEFINITION 4. Let A be a p x q matrix, and B be an m x n matrix ; then 


B = 


«o.oP 


\ flp-i.ofi 


an^-iP ^ 

a p-i,q-\B 


The p x q block matrix A Q B is called a Kronecker product . There are a number of rules for Kronecker 

products; see [25]. In this paper we will make frequent use of 


( 2 . 6 ) 


(A0B)(CPD) = (AC) (BD), 
(.4 ' B) T = A T 0 B T , 
{A;B)- [ =A- l o>B- 1 . 


Consider the following matrices, 


(2.7) 


B = B t > 0 . C = C T > 0 , D = diag(di) > 0 , 
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where B . C\ and 5® C have the structure 


" b l 


' c L 

. 0 


i 0 


, c = 


0 i 


0 

b r _ 


C R _ 


We will need the following lemmas. 

LEMMA 1. Let C and D be the M x M matrices defined in (2.7)-(2.8). If the blocks Cr.Cr have 

dimension r x r and the first and last r components in D are constant, then 

(2.9) CD = DC = C 1/2 DC l/ - > 0. 

Proof: A direct matrix multiplication leads to (2.9). □ 

Lemma 2. Let the first and last r components Dr = diag(di), k — 1 N be constants; let B have q x q 

blocks Br , Br and C have r x v blocks ( r , ( r ; see ( 2. 7)- (2.8). With .4 = B Q (' and D = diag(Dk) we have 

(2.10) AD = DA = A ]/2 DA l/2 > 0 

if the first (Di,...,D q ) and last (Dti-fa-i ). .... Dn) q -blocks in D art equal. Proof: By introducing Di = 
{I,, . ' D\) with I the identity matrix, the relations (2.6) and lemma 1, the upper left corner of 1 1) becomes 

(. B l o C)D l = B l G CD , = B L D,C = D L (B L C C) 


and 


(B l i C)D,. 




e 1 / 2 )(fi [ /2 g /v l/: ’) = (b[ /2 ( ■ r l/:; )/i,.(/t} / ' J o c 


J/-'\ 


Positive definiteness follows directly from the fact that Dr >0. □ 

3. Tlie ID difference operators. The ID difference operators that form the basis for the multidi- 
mensional difference approximations will be presented below, for more information see [17], [18]. 

3.1. The discrete differentiation operator. Let l . VU be the numerical approximations of the 
scalar quantities u and u r respectively. The approximation VU of the first derivative 


(3.1) VU = P~ l QU , P~ l Qu = T f , \T t \ = 0(Ar’ n .Ax n ) 

where |7e| = (9( Ax™ . Ax n ) means that the approximation of the differential operator is accurale to order 
in in the interior of the domain and to order n at the boundary. Typically we have n = w — 1. The 
summation-by-parts (SBP) operator PI 7 satisfies 


(3.2) 


(F.PV> = U n Vn - Wo - (VU, V) P 


where 

(3.3) (r.V)p = U T PV, P-P T . Q + Q T = D. D = diag[- 1.0 0.1] 

and 0 < < P < p„, ar ArI. Operators of the SBP type arise naturally with centered difference 

approximations; for examples see [9], [15], [12], [26]. 



In this paper we will apply the first derivative operator twice to obtain the second derivative; i.e., we 
will use 


(3.4) 


V 2 U = V(VU). u xx - P~ l QP~ l Qu = 7;, T e = 0(Ax”\Az p ) 


despite the fact that we lose two orders of magnitude in accuracy p = nr — 2 at the boundaries. The 
second derivative defined in (3.4) satisfies (3.2) with V = BVU , which is completely similar to (w, {bu x ) x ) = 
u6uj.|o — (u x ,bu x ) obtained in the continuous case. For another type of second derivative, see [17], [18]. 

3.2. The discrete integration operator. The matrix P in the SBP derivative operator is a discrete 
integration operator. 

Theorem 1. Let the difference operator V — P~ l Q defined in (3.I)-(3.3) exists on the interval -1 < 
z < T There the matrix P is an integration operator which satisfies 


(3.5) 


J (uv) T dx = U T PVV + ( VU) T PV 


where U, V are the projections of the continuous functions u, v onto the grid. Proof: U T PVV + (VV) T PV = 

(cv)IS f = («w)lli = /i 1 («»M*. □ 

It is also possible to prove the following theorem. 

Theorem 2. Let the difference operator P~ X Q defined in (3. l)-(3.3) exist on the interval —1 < x < 1 
and be accurate to order m. Then . the matrix P is an integration operator which satisfies 


(3.6) 


L 


uv dx ~ U T PV + 


')■ 


where U,V are the projections of the continuous functions u,v onto the grid. Sketch of Proof: The proof 
has two parts. The first part shows that (3.6) holds for general polynomials. Next, Weistrass interpolation 
theorem (see [28]) is used to show that it also holds for continuous functions. 

4. The continuous problem. The definition and specification of the continuous problem is done 
with Cartesian coordinates. After the transformation to curvilinear coordinates we check that the essential 
mathematical properties are preserved. 

4.1. Cartesian coordinates. The two-dimensional (2D) linear problem considered in this paper is 


(4.1) 


Uf -f- P X + Lly = /?, [j 1 ,?/] G / > 0 

« = /, [;r,t/]e£T t = 0 

Lit ~ //, [z, y] G t > 0. 

where h,f,g are the data of the problem, L is the boundary operator, and 

F = F / + F V ', F f =am, F v = -(b u u x + b l2 u. y )< 

Ct = Ct 1 -f G l , = fli'U, f/' — — ( b‘2 1 U r -j- b r 22lly). 

The coefficients a-iMj are known functions of z, y y and t. For simplicity we have chosen fi = [x, y] E 
[—1, 1] x [0, 1]; see figure 4.1. For future reference we also introduce 

(4.3) a = («!,«*_>). F = (i\G). n — (ni,» 2 ), 

where n is the outward pointing unit normal on SQ. 


(4.2) 


Ecpiation (4.1) can be thought- of as a model for the Euler, Navier- Stokes, or Maxwell's equations. 
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FlG. 4.1. The computational domain. 


4.1.1. The energy estimate. Well- posed ness of (4.1) requires that we ran obtain an energy estimate. 

(u,v)=: j j uv dxdy, (u, it) = ||u||~\ 

(4-5) (u,v)EM m = / (**'') E,w dy, (ti, u) E W - ||i/||^ n 

./(> 


Let 

(4.4) 


(4.6) (u,v)n,s = j (uv) m \\sdx, |MLy i5 = (u. u)^ s ||w||iv, 

denote the Lo scalar product, the Lo norm, the boundary scalar products, and the boundary norms respec- 
tively. The subscripts E, IT, A\ S refer to the EAST, WEST, NORTH, and SOUTH boundaries (see figure 
4.1). 

The energy method applied to (4.1) leads to 


IMIr = - [(«. F 1 + 2 F v )e - (U. F 1 + 2F V )„ ] 

EAST-WEST 

- [(«, G 1 + 2 G v ) n - (u, G ] + 2G V ).s] 

NORTH-SOUTH 

- [(«. F‘) - («„ F 1 ) + (u.G 1 ) - (uy,G')] + [(«, h) + (h, «)] 

S v / N v ' 

GROWTH 1 GROWTH 2 

(4.7) - [(«„ F v ) + (F v .u r ) + (u v .G v ) + (G v . u v )} . 

s v " 

DISSIPATION 

GROWTH 1 (GR.l) and GROWTH 2 (GR2) in (4.7) will lead to a growth or decay in ||u|| 2 , but will not 
affect well-posedness. Note that for constant coefficient problems GR.l is zero. To bound ||u||* in time, the 
first two terms must be bounded using the correct boundary conditions and the DISSIPATION (DI) must 
have I In' right sign. 

4.1.2. The dissipation. For a correct sign of Dl (see (4.2) and (4.7)) the eigenvalues of B T B J where 



b\\ > 0 , boo > 0 , ^ 11^22 — [ ] > 0 


must be positive; therefore, 

(4.9) 

is a requirement for well-posed ness. 

4.1.3. Boundary conditions. By integrating the “cross derivative” components of the first two terms 
in (4.7) one obtains 

- f o u(F ! + 2F v )\t\dy- J ^ u{G* + 2G v )\ l 0 dx = (6 12 + 6ai)« a |±ilo 

- J w((ai + (6i 2 ) y )« - - j u((a-> + (&2i)x)« - 2 bn- 2 U y )\odx . 

The term (612 + ^ 2 i)« 2 |lilo involves the point values in the four corners of the computational domain. These 
point values cannot be estimated (unless they are artificially specified). Boundary conditions where the 
viscous fluxes (F' ,G l ) are specified avoid that difficulty and lead to an energy estimate; therefore one 
should specify the total flux (F = F 1 4 - F v , G = G 1 + G v ) or the viscous flux at the boundaries. 

Consider the first two terms in (4.7) and recall the definition (4.3). At x = ±1 with n = (±1,0) we have 
the boundary conditions 


(4.10) 


while 


-a i(- 

- 1 . (/,<) 

= a 

■ n < 0 

F 

* 71 

= F 

= /'H 

(y-t), 

— «i(" 


— a 

• 77 > 0 

F v 

• ?7 

= F v 

— F v 

— r w 


+«i(+i, y< 0 

= a 

• n > 0 

F v 

• ?7 

= F v 

— /A' 

“ r A 

(iM). 

+«i (+1 , y, t) 

= a 

• n < 0 

F 

■ i7 

= F 

= Fs(y,t), 

—ao(. 

r, 0, t) 

= a • 

n < 0 

F • 

?7 

= C 

= Gs{ 

•M). 

— do ( . 

r.O J) 

— a • 

71 > 0 

F v ■ 

?7 

= G v 

= GU 

x,t). 

+a 2 (. 

r,U) 

= d * 

/7 > 0 

F v • 

?7 

= G v 

= <?&( 

Xj), 

+a 2 (. 

r.l.f) 

= a • 

n < 0 

F- 

/7 

= r; 

- (F\ ( 

•M), 


(4.11) 


are used at y = 0, 1 where n = (0,Tl) . The boundary conditions (4.10), (4.11) can also be formulated in 
the following more general way. 


(4.12) 


a * n < 0 => F ■ n — Fsn • n 

a n > 0 =$ 


?v . n - 


sn 


Boundary conditions of the type (4.10), (4.11), (4.12) have been derived in [29] and [ 20 ] for the Navier-Stokes 
equations. 

Let us consider the EAST boundary in detail, and let us assume that a\ is positive at y = 0, becomes 
negative at y = y 0 and remains negative until y — 1 . Inserting the boundary conditions (4.10) into (4.7) 
yields 

r l , /‘.Vo 

- j ii(F ! +2F V )\£(ly = - |«i ( 1 , j/, 0 1« 2 + 2 uF E dy 

Jo ./(I 

ry° 

+ / (1, y, t) ir — 2uFt'dy = — / ja 1 |m“ + 2 «Fj ; dy 

J y Q J 0 

- I |«i|«" + '2uFEdy = - I |«i|ir + 2 </(rr 1 F l E + F E )dy , 

J y u J 0 


(4.13) 



where cr\ — (1 — |fli|/cii)/2. The requirement that a \ changes sign in the manner described above can be 
relaxed, so (4.13) is a generally valid formula. An entirely similar procedure at the other (WEST, NORTH, 
SOUTH) boundaries yields the final result: 

-(«. F' + 2 F v )e = -(«, |«i|w)e - 2(u, Fe)b- 
(«. F 1 + 2 F v )w = -(«, M«)i v - 2(«. F„ )„ , 

~(u,G' + 26' 1 ).\- = — (m, |« 2 I m ) at - 2 (u.Gn)x. 

(4.14) ( u.G 1 + 2 O l ),s - —(it, |a 2 |w).s — ‘2(u,Gs)s- 

where 

FE=<r l F B + [i-<ri)F%, <n = (1 - |«i|/a t )/2, 

F\\ = PzFw + (1 — 1 T 3 ) = — (1 + | a 1 1 / o 1 ) / 2 , 

Gn — <r 5 Gs + (1 - <t 5 )G v n . <t 5 = (1 — |a 2 |/a 2 )/2, 

(4.15) Gs = <t 7 6's + ( 1 - <t 7 )G'' s , <t 7 = -( 1 + |o 2 |/a 2 )/2. 


(4.16) 


where 


Inserting the relations (4.14) and (4.15) into (4.7) leads to 

1 


w|lt< S — IIC/II/ + CR1 + GR2+DI, 
r=E,w,N,s T ' 1 


r]E = 


in K 

In " 2d y 


mv 


,/u |ai|« 2 rfy , 
Jo u-dy 


\x-~ 1 1 


fl, lazlirdx fi, |a 2 |«"£/j.‘ 

V.v = — n — — ly=i ’ = — 71 — — lw= n • 

J_ l u-dx J_ j u-dx 

The parameters ije , ?/h'i 7 /s are strictly positive if a^ao are zero for a finite number of points. 

Time-integration of (4.16) leads to an energy estimate of the form (2.3) if (4.9) holds. Provided that 
a solution exists (can be shown by using the Laplace- transform technique or via difference approximations; 
see [30] and [31]), we can conclude that the following theorem holds. 

THEOREM 3. Problem (4 A), (4-10), (4 At) is strongly well posed. 

4.2. Curvilinear coordinates. In this Section we consider problem (4.1) on a curvilinear domain. By 
introducing the transformation t — r, x = .r(£, ?/, r), y — ?/(£, ?;. r) and it's inverse 


(4-17) 


T = /, 


i(*\yA), y - t}(x<yA), 


we obtain the transformed equation 

(4.18) u )t A (J (m« A ZjtF A £ y G))z A (J(yt 11 + i]xF A %C)) r/ — Jh A RR 

where 


IUIS = u(J t + (J£tU + (J*hU) + F(VirU + A G((JZyU A 



Using the metric relations causes the term RHS to vanish: 


Jtt 

= (x t) y T - x T y t) ) 

•/Vr 

= (ys*T ■ 

- x (.Vt) 

d s.r 

= 

J£y 

= -x„ 


J*)r 

= -«e 

J>ly 

= 


J 


J 

= (6r>Jy - 

~ £. y 1 h ) 


In this paper we will consider the steady version of (4.17), i.e., £ = £{x,y),ii = r)(x,y). The new 
transformed problem becomes 

Ju T + (F)f -f ((7) r j = /?, [£, 7}\ G t > 0, 

(4.20) u = /, [f , Tj] G ft, . r = 0, 

Lu = g, [£, 7]] G <Jft, r > 0 T 

where h = Jh, f,g are the data of the problem and ft = [£, r;] G [— 1, 1] x [0, 1]. The new' transformed fluxes 
are 


F = J(F • V£) = F 1 + F v , F^aiu, F v = -[6„« { + 6 12 u„], 
(4.21) G=J(FVi]) = G , +G v , G , = a 2 u, G v = -[b 21 u ( + 6 22 u„], 


where 

(4.22) 


a x = Ja ■ 6 n = J Vi ■ BVt, 6, 2 = ■ BVy 

a 2 = Ja ■ V>/ & 2 i = JVi) ■ _BV£ 6 22 = JV// • BVrj 


and B is given in (4.8). 

4.2.1. The energy method. Let 


(4.23) 


(4.24) 


(w.v)j 


U> 


Jdtdr /, ||w||j = (w, u)j, 


(«•<')=/ J («*■’) |M| 2 = (u,u), 


(4.25) 


(«, ?')e.iv 



IMIe.jv 


(u> «)k.h • 


(4.26) 


(«, t’)tV,5 = 





ll«ll?V,5 = («.«)tV.S 


denote the weighted L 2 scalar product and norm, the L-> scalar product and norm, the boundary scalar 
products, and boundary norms respectively. The subscripts E, IV, A 7 , 5 refer to the EAST, WEST, NORTH 
and SOUTH boundaries as in figure 4.1, with x,y replaced by V //. 

The equation corresponding to (4.7) becomes, 


(Il«ll5)r = - [(». F^ + 2F v )e — (n. V f + 2F 1 )„ ] 
EAST-WEST 

— [(«., Cr^ -f- 2(7 ' ),V — ( “H 2(r' ) S'] 

' 

NORTH-SOUTH 


8 


(4.27) 


- [(«, F{) ~ («{ , F‘ f ) + (u. G' n ) - («„, G 1 )] + [(», h) + (/i,ii)] 

GROWTH 1 GROWTH 2 

- [(«(. F v ) + (F v .ti f ) + K-G'j + (G r , «,)] . 


DISSIPATION 


Precisely as in the Cartesian case, GR1 and GR2 in (4.27) can lead to a growth or decay in \\u\\j, blit will 
not affect well- posed ness. The metric relations (4.19) show that in the curvilinear case too, GR1 vanishes 
for constant coefficient problems. Just as in the Cartesian case, we need to assure that DI has the right sign. 

4.2.2. The dissipation. Similar conditions as in the ( 'artesian case for positive eigenvalues also apply 
in the curvilinear case; see (4.27). Thus we must show that 


(4.28) 


b }\ > 0, ^22 ^ O' 611602 — 




to assure the right sign on DI. The conditions (4.28) hold, since (4.22), (4.8), and (4.9) lead to 26u = 
JV£ T {B + B t )V£ > 0,26o2 = JVt) T (B + B T )Vi) > 0, and 


61 1600 — 


612 -f 6*21 
_ 


= 6] 1600 


6 jo *T 60 


> 0. 


4.2.3. Boundary conditions. Consider the first two terms in (4.27) and recall the definitions (4.3). 
The outward pointing unit normal on is. 


(4.29) 


»(£ = ±1 .»?) = = 0, 1) - 


iv«r 


|V»/| ’ 


where |V£| = ^Zr + and |V»/| = ^Ji)]. + f/jj. The boundary conditions leading to an energy esti 


stimate 


become 


(4.30) 


-oi(-l.t/,<) = —Ja ■ < 0 JFVZ = 


— «i ( — 1 , ?/. t) = 
+ai (+!,>;./) = 


-Ja V£>0 JF V -VZ 
JaVZ> 0 JF V V£ 


F = l-\ v ( V J), 

F y = r? y ( v j), 
F v = FjAq.l). 


at £ = ±1, while 


(4.31) 


+«1 (+],?/,/) = 

Ja • V£ < 0 

JFVZ = F — 

Fe(’iA) 

-«2 (COT) = 

— JaVif< 0 

JFVr ) = G = 

Gs(Ft)- 

-r/ L ,(C0J) = 

— Ja V i} > 0 

JF V Vi) = G v = 

G'pCO, 

+a 2 (ZA.i) = 

Ja Vi) > 0 

JF V ■ Vq = G v = 

G) V (C/). 

+a?(ZA.t) = 

Ja • Vi) < 0 

JFVq = G = 

G'vfCO 


should be used at // = 0. 1. The boundary conditions (4.30). (4.31) can also be formulated as in (4.12). 
The same procedure as in the (’artesian case leads to the estimate 


(4.32) 

where 


m||5) t < V — 1|/>||7 +GR1 +GR2+DI. 

• * 1 1 r 


7=£\M \ ,.S 


»// 


Ft: = & 1 Ft: + ( I - )F X E . tr { = -(1 - |«i |/«i ) 
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(4.33) 

and 


F\v = + (1 - <t 3 )F,V, <t 3 = ~ ( 1 + |«i|/fi|), 

Gn = (T $ (J tv + (1 — (T 5 ) G , <r 5 = -(1 — Idoj/ao), 

C*s = <77^/5 4- (1 — (77)6x5, 0-7 = — -(1 -f |a 2 |/«2)i 


So l«il« 2rf, ?i 
= f l ■>, U= 


So u ~ d n 


So l« 

W = — 


| vrdr) 


So u ' dr l 


k=- 


_ S-l I«2|m 2 ^ 
/_r« 2 ^ 


S- 1 l a 3l« 2 ^ 
S-i u ~ d Z 


1*7 = 0- 


The parameters ?/n , On, f)s are strictly positive if d \ , a 2 are zero for a finite number of points. 

Time-integration of the estimate (4.32) leads to an energy estimate of the form (2.3) if (4.28) holds. 
Provided that a solution exists we can conclude that the following theorem holds. 

Theorem 4. Problem (4.20), ( 4 . SO), (4-SI) is strongly well posed . 


4.3. Interface conditions. Boundary and interface conditions of the flux type (see (4.10), (4.11), 
(4.30), and (4.31)) require extra careful treatment; see [27] for an example. 

4.3.1. Interface conditions in the curvilinear case. To apply the SAT technique [16] on the fluxes 
at an interface between two blocks with different coordinate transformations and matching gridlines (see [17], 
[18] for the one-dimensional treatment) requires that we identify the continuous part. Matching gridlines at 
£ = £0 = const implies 


( 4 - 34 ) (*?)l # (^) 2 , (^7,)i = (%)i = [yr,h 

while we have 

U-35) = (*e)2. (y*h = (%) 2. (^v)i 7 L (xvh, (y n )i ^ {y n ) 2 

at T) = 7/o — const . The subscripts 1,2 refer to the two coordinate transformations. 
Equations (4.21), (4.19), and (4.34), (4.35) immediately lead to the conclusion that 

( 4 ' 3G ) fiKo, T ) = ^ 2(61 T )< Gi(£o, t) ± G 2 (£ 0 . r), 


( 4: * 7 ) Fi(S<*IQ'T) 7 ^ F'ACO OiT), Cn{£, i)q,t) - G 2(£.r/ 0 ,r); 

i.e., F is continuous across £ = ccwsf while G is continuous across 1 / =r const . 

4.3.2. Interface conditions and vanishing wave speeds. Another problem with flux-interface con- 
ditions appears when the wave speed a goes to zero. Consider the two constant coefficient problems 

u t + F(u) r = 0, —L < x < 0 and v t -f F(r) x = 0, 0 < x < L, 

where F(w) — rnr + F' (tr), F v (w) — — otc x . Both problems have homogeneous outer boundary conditions 
at |.r | = L and zero initial data, and they are connected through interface conditions at x = 0. We 


JO 


will compare the effects of flux-interface conditions (F(u) = F(v),F x («) = F x (c)) and variable-interface 
conditions (u = r, u x = c r ) on the solutions. 

By transforming the problem for v on [0, +L] onto [— L, 0] via the transformation x — £, and then 
replacing £ with i\ we obtain 


(4.38) 


+ A V’x 


t > 0, 

-L < x < 0, 

%■ 

= 0, 

1 = 0, 

- L < j- < 0, 

B-lH’ 

= o, 

i > 0, 

.r = -L, 

Bat 

= 0, 

/ > 0, 

X = 0, 


where v — (u,t>) T ,A = diag(a,—a), and B-if = 0 denotes the outer boundary conditions. Buy' — 0 
represents the transformed interface conditions 


(4.39) an — eu. x — av + ay, — e u x — +(v x or u — r, f/j. = — r r . 

We will treat (4.38) as a half-plane problem, which means that we let L — > oc and replace the influence of 
B-l by only admitting bounded solutions as x — > — oc. 

The Laplace-transform technique applied to (4.38) leads to 

u(x , s) = <T\ (s) exp (Ki(ti)x), v(x<s) — <r 2 (s) exp (k 2 (s)c) 


where s is the dual variable with respect to time and 


a I a s 

Kl - + Tt + V ( Tc y + 7' 


a 




s 

( 


Note that both u, and v decay away from the boundary x ~ 0. 

The interface conditions (4.39) lead to the equation E(s)d — 0 where ff — (^i,^)^- A well -posed 
bounded solution is obtained only if det(E(s)) ^ 0 for :)?(s) > 0. The flux-interface conditions in (4.39 lead 
to 


=> det(E(s)) = -2 ( a ] J{Ey'+^ 

dft(E(s))= + 

Obviously the flux-interface conditions can lead to unbounded growth for vanishing wave speeds, because 
det(E) a _ ¥0 = 0 independent of s. The variable-interface conditions, on the other hand, lead to a well-posed 
problem since dr t( E) a _^ {) = 2\/(s/(). 

A similar analysis of the flux-boundary condition au — (u x = 0 for the single domain yields drt(E(s)) = 
a/2 -f \J(a/ 2)~ + s ( . Consequently, the problem with unbounded growt h for vanishing wave speed does not 
exist in the boundary condition case because dft(E) a _ ¥{] = \J ( sr ) . 

5. The discrete problem on a curvilinear mesh. In the rest of this paper we will consider the 
transformed problem (4.20), (4.30). (4.31). Note that problem (4.1), (4. 10), (4. 1 1 ) corresponds to the special 
case where r = / ,£ = x and ?/ = y. For notational simplicity, we ignore the "hat" notation for the fluxes and 
transformed coefficients introduced in (4.2U)-(4.22). 


/ a — ck\ — a — €K r > 

(4.40) E(s) = 

V — CK 1 — €Ko 

while the variable-interface conditions leads to 

(4.41) E(s)=l 1 _1 ) 

V «1 «2 / 


n 



NORTH 



£ = -l * = +l 

SOUTH 

Fig. 5.1. The single- domain case in transformed space. 


Let. the N x N matrix and the M x M matrix P r , be the ID positive definite matrices defined in 
Section 3.1. A product av is arranged discretely (where av ^ .4 V ) in the following way (see Figure 5.1): 



respectively. The subscripts E, W, N, S refers to the EAST, WEST, NORTH, and SOUTH boundaries (see 
figure 5.1). At this point we also define the restriction of U to the boundaries: 

(5-4) Fejv = {Je w 0 Aj)Ui Uv,s = ( A - Jx.sW- 

b.l. The norms in the transformed, problem. The norms and scalar-products corresponding to 
(4.23)-(4.26) are 


5.5) 

II 

' — i 


'• (r.r 

h = IIHI: 

o.fi) 

(C V) = 

n>u. 

■ (C.U) 

=11^. 

5-7) 

(U,V) E ,w = U t (Jb.w P v )V = 

■ 

u„u E .„ , 

\m\hv 
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SOUTH 

FlG. 5.2. Schematic showing J requirements necessary for A iJ to be a norm: r — 4 and q = 2. 

.(5.8) (U, U),v,5 = U T (P< © Jn,s)V = l|f'|| 2 v,5 = (E l T ).\ i.s- 

Obviously, the relations (5-6)- (5-8) define norms since and P 7) are positive definite matrices. What about 
(P* P*)J in (5.5)? 

The metric scalar J is defined in (4.19). In matrix formulation we have 

(5.9) J — diag(Ji) y i = 1 . N J,- = diag(Jij)J = 1 A/. 

We need the following theorem. 

Theorem 5. Let M = ■ P^. If the first and last r components in Jj are constants and the first 

(Ju--..Jq) and last Jn) q blocks in J are equal, then MJ is a norm. Proof: Lemma 2, with 

B — P^C = Prj, and I) — J . with J defined in (5.9), leads directly to MJ — JM — M [ ^ J JM [ ^ 2 > 0 □ 

The requirements in theorem 5 are illustrated in figure 5.2 (for r = 4 and q — 2) where rxq values of J are 
equal in the corners and r, q values of J are constant normal to an rj — const ^ = const boundary respectively. 
A curvilinear transformation with such a J close to £12 is called volume preserving and guarantees that MJ 
is a norm. 

Remark. The conditions in theorem 5 (i.e., that J must be constant in the first r/, r points normal and 
adjacent to the boundary £12) can be t hought of as theoretically ideal conditions. In practise one approaches 
the ideal condition with increasing resolution on a smooth mesh close to the boundary because 

J(iJ) ~ J( OJ) = 4(0, f]j)(iAO + O(Ae), i = 1 q. 

J(iJ) - J(i % ()) = J,(6,U)(jAt;) +0(Aiy 2 ), j = 1 r. 

where it is assumed that J(0. j), «/(■/. 0) are the values of J at the boundaries. This process is illustrated in 
figure 5.3, where the minimum eigenvalue of PD+ DP as a function of increasing resolution is shown. The 
minimum eigenvalue goes from a negative value for large A.r to a positive one for small Aj\ 

5.2. The single-domain problem. The discrete formulation of (4.20), (4.30). (4.31) with the SAT 
technique [1(5] for incorporating flux boundary conditions is 

(5.10) JV T + lh F + A, = h + B( \ P(0) = /. 

where the continuous derivatives are approximated with 

(5.11) d^f = (p-'q v i„)F, D„(; = (i r r,;'Q n )C 
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GRID REFINEMENT ON (PD +DP) 



D n = 1 + 100 sin(alpha x,) 2 
P = 6th-order norm 


negative eigenvalues 

positive eigenvalues 


J 1 1 i i i i i » 

100 200 

Gridpoints 


FlG. 5.3. Minimum eigenvalue of PD + DP as a function of Ax. 


and 


BC = (Pp J E ® 7,Ei)(F - F e ) + (PpJ E © /„E 2 )(F l - F%) 

+ (PpJw © /,E 3 )(F - F w ) + ( PpJw 0 /,S 4 )(F v ' - F,\) 

+ ( ©5 © P n 1 Jn)(G — GaO + ( I ^ ©f; ® Pp Jj\)(G x — Ci ^r) 

(5.12) + (/ c Et © F“ l J.s )(6’ - G s ) + (/ € S« © P" 1 Js)(G v - <?£). 

The A r x A r matrix and the M x M matrix Q n are defined in Section 3.1. Fluxes with subscripts E , IF, A r , S’ 
are boundary data. The matrices Si — E 8 will be determined below. 

5.2.1. The energy method. Multiplying (5.10) from the left with V T (P$ G P n ), introducing the 
notation M = P$ O Pr), and adding the transpose of the equation leads to 


)r = -[U 1 (Qt < P„)F + F 1 (Qj e P tl )U] + [U T Mh 4- F r Ml 


(5.13) 


A GR2 

+ -W T (Pt : Qr,)G + a T (P ( (o <^)P]+sr + (BT) T 

s. ^ 

B 


where 


FF = o P r ,Si)(F - Fe) + U T (J E 0 P„S,)(F r - f£) 
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+ U T (J W O P„Z 3 )(F - F w ) + U T (J W O P^4)(F v - F,V) 

+ r rT (P«S 5 O Jjv)(G - 6 a-) + F t (P<S 6 0 Jn)(G v - C£) 

(5.14) + U t (P<Z 7 r. Js)(G - G„) + C T (P £ S 8 O' J S )(G V - G'V ). 

In (5.13) we have assumed that the metric transformation is such that MJ is a norm. Note that the flux 
terms with subscripts are given data. 

The notations and abbreviations 

(5.15) Qi + Qj = B ( = J e - Jw, Q n + Q% = B ’i = j n ~ J *' 
will be used to expand the A and B in (5.13). We obtain 

A+B = - [U t (B ( C P n )(F ! + 2 F v ) + (F 1 + '2F V ) T (B ( O P v )l >]/ 2 

V ■ - -V 

E-W 

- [F t (Pz O Br,)(G' +26' ) + (O' + 2 G V ) T [P ( O B v )U]/2 

s / 

N-S 

- {[(U. D ( F l ) + (D^F 1 , (' )] - [{F 1 , D<U) + (D £ C. F')]}/2 

S J - s 

GR1 

- {[(V, DrfG 1 ) + (D V G'.U)] - [(6 7 . D„U) + (D n U, 6 / )]}/2 

' ^ ' 

GR1 

(5.16) + [(D ( U, F v ) + (F l \ D ( U) + (D„U),G V ) + (6 r . D„U)] . 

\ ^ - 

DI 

Note the close similarity of the discrete energy estimate (5.13). (5.14), (5.1(j) with the corresponding 
continuous one; see (4.27). Just as in the continuous case GR1 and GR2 will at most create an exponential 
time growth. To obtain an energy estimate we must determine under what conditions the dissipation (DI) is 
negative definite and which values we must assign to the matrices Xb — H s to obtain bounded contributions 
from the boundary. 


5.2.2. The numerical dissipation. The DI is 

(5.17) DI = (D ( U) t AIF v + ( F v ) r .UDzU + (D„F) T MG y + (G v ) T MD v r. 


The relationship between the gradients and the fluxes in the continuous case are given in (4.21). In matrix 
formulation for the discrete case we have 

(5.18) 


’ F v ' 


’ B n 

B \ 2 


' /V ' 

G v 


_ B, i 

B 22 


D„F 


where Bki(Lj) = //;). By introducing (5.18) into (5.17) we get 

B\\M + M B\ \ B 2] M 4- MBn 
B V1 M + M B 2 1 B T2 M + MB T1 


(5.19) 


DI — — 


' D<r ' 


D„C 




-—1 


D„r J 


At this point we need to define M = P £ F y J} in greater detail: we have 


• 


■ 

1 

0 


' "n 

1 

0 

(5.20) 


0 

I 

"f . 

• p„ = 


0 

1 

//* 

V J 
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We will need the following lemma. 

Lemma 3. If the blocks have the size q x q } the blocks have the si ze r x r, and the 

matrices B^i in (5.18) are constant in the first q,r points normal and adjacent to the boundary SSI. then the 
dissipation DI defined in (5.19) is negative definite. Proof: Lemma 2 leads to 

01 1 A/ T A I B\\ -02i A/ -j- M 0i2 

By>AI -+■ M 021 022 A'l + M 022 

M i/2 o 1 r 2 B n B 1S + B 21 1 [ M 1 ' 2 0 1 

0 A/ 1 / 2 J [ B 12 + B 21 2 B 22 0 A/ 1 / 2 >0 Q 

Remark. The conditions in lemma 3 (i.e., that the matrices B^i in (5.18) are constant- in the first q,r 
points normal and adjacent to the boundary <5f}) can be thought of as theoretically ideal conditions. In 
practise one approaches the ideal condition with increasing resolution, smooth coefficients and a smooth 
mesh; see the Remark on J in Section 5T. 

5.2.3. Stability. To obtain an energy estimate (given the negative contribution of the DI on the right- 
hand-side of (5.16)) we must assign values to the matrices Ei - E 8 in order to obtain a bounded boundary 
contribution. 

Let us start by estimating the terms at the EAST boundary. We have 

BT e = - {U T [P„(I/2 - SOJF' + (F'f[(I/ 2 - Ef )P„]F) 

- {U T [P n (I - Ei - E 2 )]F*' + (F V ') T [(/ - Sf - Ej)P,]C} 

(5-21) -[ U T p ti F E + (F E ) T P n Ul 

where F E = E t F E + • 

Obviously, the terms involving the viscous fluxes must be removed. This yields = / - S x . By 
observing that F 1 = A E U where \ E = diag((ai) Nj ) (see (4.21) for a definition of a,) we obtain, 

BT e = - U T [Pn(I/ 2 - E x )]A e + A £ (//2 - Sf )P„]U 
(5-22) - [U r P„F E + (F E ) T P v V]. 

Now we choose Ei such that (1/2 — EiJAf: = |A^|/2. This choice and an entirely similar. procedure at the 
other boundaries yields 

(5-23) (||C||3) r < —11^/11/ + OR1 + GR2+DI, 

I=E,\\\N,S 

where 

F e = Ei F e -f (I y - Ei )F x e , Ei = (I y - |A^|A“ 1 )/2, 

F W = S 3 F» + (Iy -r 3 )/,V. S 3 = -(Iy + |An- jA n -‘ )/2, 

6',v = SgGjv + (I T - ^5 )G\’. s 5 = (/.,. - |A ; v|A j y 1 )/2, 

(5-24) G s = S 7 C.s- + (I r - S 7 )G.V, E 7 = -(/, + |A 5 |A" 1 )/2, 

(5.25) 
and 


— I —V 1 V. — _/ V' V I V v j v 

— Iy ~j\ ♦ -^4 — iy — — Ij “ —'5 i -<8 — — Ar — E7, 


l(C|A,|r), + (|A,|CC),] , = flr¥ , 

1 2 (r : , F)z ' ■ , ’ 5- 
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The similarity of the discrete energy estimate (5.23) with the corresponding continuous one (see (4.16) 
and (4.32)) implies strict stability. Time-integration of (5.23) leads to an estimate of the form (2.4) if the 
DI has the right sign, i.e., if lemma 3 holds. We can conclude that the following theorem holds. 

Theorem 6. The approximation (5. 10) of the problem (4-20), (4-20), (4-21) is both strictly and strongly 
stable if lemma 3 holds and Hi — Hh are given by (5.24) (in ^ (5.25). 

5.3. The ID multiple domain problem revisited. Before we consider the 2D multiple-domain 
problem, let us once more look at the ID multiple-domain problem considered in [17], [18]. 

5.3.1. Derivation of the Q-formulation for interface problems. Consider the following hyperbolic 
interface problem 

(5.26) lit + u x = 0, — 1 < x < 0, and v { + v r =0, 0 < x < 1 

augmented with suitable initial and boundary data and the interface condition u = v at x — 0. The 

straight forward approximation of (5.26) is 

Ut + PpQ L l J = erVne.v - V'o)r.v ) 

(5.27) V t + Pr'QrV = — Us)e 0 ) 

where V = (U a U N ) T .r N = (0 0, l) T , V = (I n, .... V M ) r .e 0 = (1, 0...,0) T . 

The boundary terms from the left (L) and right (/?) outer boundaries are ignored. The formulation 

(5.27) can also be written in the following way: 



The 2 x 2 blocks of Q sk and D corresponding to the nonzero elements in H are 

0 d - - 1 “ ' 2<j l vl + ^r 

-(<t l-ctr) 0 1 2 [ (tl + —\ — *2<jk 

In the sequel, the “tilde" sign will indicate the 4 x 4 block that couples the solutions in the left and right 
domains. Equation (5.28) now becomes 



(5.30) P\Y, + (Q' k + D)W = 0. 

In [17] it was shown that (5.27) is conservative if an = (Tl — J . By introducing this condition in Q' k and D 
we obtain the final form of the difference operator 


(5.:u) 




where a ~ 1/2 — <t l . 

The formulation (5.30), (5.31) hereafter referred to as the Q-formulation is a rearranged form of the orig- 
inal formulation (5.27). However, the Q-formulation simplifies and even extends the possibility to formulate 
suitable penalty terms for second order derivatives. 

5.3.2. The Q-formulation for advection-diffusion interface problems. Consider 


(5-32) u t + F(u) x =0, “1 < •*- < 0, and v t -f F(t?) r = 0, 0 < x < 1 


where F(w) — a(x,t)w — ew T augmented with suitable initial, boundary, and interface conditions. An 
approximation of (5.32) using the Q-formulation is 


(5.33) 


PW t + (Q* k )(AW) - e(Q sk + D 2 )P~ l (Q sk + D 3 ) = D\ W 


where IT — (U,V) T and P = diag(Pi , Pr). The matrix A has the values of a(xi J) on the diagonal. The 
operator Q' k is defined in the previous Section, and 


(5.34) 


D t = <Ji 


1 -1 

-1 1 


* = 1,2,3 


as in (5.31). The dissipation D\ is formulated as acting on TT , which is a more general formulation that 
includes penalty on the flux (< 7 ! = <ra(0,/)) as well as penalty on the variables. 

We can now prove 

Theorem 7. The approximation (5.33), (5.34) °f the problem (5.32) with the choices 


(5.35) 


<T\ <0, <72—0, <7 3 — 0 


is conservative and stable . 

Pioof: The energy method applied on (5.33) leads to 


II W'llf = (VW, AW) - (V(AW), W) - 2 e(VW, VW) - W T B(AW - 2e VW) +IT 

S ^ V V ^ ^ V. y 

GRl Dl BT 

where VW — P~ l Q sk W and the interface terms IT are defined as 


(5.36) 


w " 

T 

2Di +2 cD-,P~ 1 D 3 

({D‘2 — D3) 


w 

VW 

0 

((Dn — D 3 ) 

0 


VW 


o 


The growth (GRl), the dissipation (Dl) and the ordinary boundary terms (BT) match the terms in the 
corresponding continuous estimate perfectly. The choices (5.35) makes the term IT maximally negative 
definite and leads to stability. The approximation (5.33) can now can be written 


(5.37) 


PW t + Q* k (AW - eP~ l Q $k W) = D { IT, 


which leads to conservation. □ 

It remains to identify the penalty terms and corresponding interface conditions and find out whether 
(5.37) is sufficiently accurate. The equation 

PW, + Q(AW - iP-'QW) = (Di + (Q - + f (Q sk P~ l Q' k - QP~ l Q)W 
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FlG. 5.4. The multiple domain mesh in transformed space . 

is a formulation of (5.37) in the usual penalty form. Obviously, the penalty terms (denoted by PT) indicate 
that the one-sided first and second derivatives are replaced by first and second derivatives involving the 
adjacent domain, and that dissipation is introduced via D\. 

By introducing 


> 

II 

1 

1 -1 * 

, A = 

-{2ai + a) 

0 

2 

1 -1 


0 

(2<T! - «) 


we get PT = AAM +< ( A(P [ QW) -f Q 5 ^AIU); this result shows that the corresponding interface conditions 
are u = v and u x — v T . Also, because 

AcS> = (0,0) T , A(P~ 1 Q</>) = (0( A* 1 , A*Jf 1 ) , 0{Ax r L ~ 1 . Ax r ^ 1 ) ) T 

where o is a smooth function and r is the order of the approximation at inner points, the approximation 

(5.37) is accurate enough. 

5.4. The 2D multiple-domain problem. In this Section, an interface at £ ~ 0 with matching 
gridlines (see Figure 5.4) is considered. Matching gridlines implies that the number of points in the ;/ 
direction (M) is the same on both sides of £ = 0. We will also assume that = P ^ — P v . which implies 
that Qfj = — Q T} . Note that, in general, the difference operators , D ^ can be different in the left and 

right domains and that A£^ ^ A£/? and Nl ^ A/?. 

A multiple-domain Q-formulation of the problem (4.20), (4.30), (4.31) is 

(5.38) JVY, + D\ k F + D tl G = M~ 1 (D SP„)W* + h + BC, U(0) = / 

where U = (/’, \ ' ) T . The solutions in the left ( L) and right (/?) domains are denoted respectively bv V and 
W and 

(5.39) Df = U-'iQ? P, } ). D n = M~'(Pi . Q„). 

In (5.38), BO denotes the boundary conditions in (5.10) at the NORTH, FAST, SOUTH. WEST boundaries 
in penalty form, h is the forcing function, / the initial data, and F — F 1 + O' .(7 = ( j 1 + (V' the fluxes 
where 

(5,10) i.coF 1 T,U. F v = -(R u Dfir+ H v ,D„\V), 
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(5.41) G‘ = A 2 W, G v = -(BniDfW + B 23 D V W). 

The remaining definitions and notations use.d in (5.38) are Q^ k = + A, 


(5.42) 


(5.43) 


M = 


A = 


M l 0 
0 Mr 


J = 


0 0 

A 

0 0 


Jl o 
0 Jr 


D = 


■ Q ( = 

0 
0 


D 


Qf 0 

o Qf 

0 
0 


p t 

0 1 

1 

’ 1 

-1 ' 


1 

-1 " 

0 

1 

’ A = -2 

1 

-1 

, D = 

-1 

1 


.(5.44) 


The matrix coefficient E will be determined by stability requirements. 


5.4.1. Conservation. The Q-formulation automatically leads to conservation: 

Theorem 8. The approximation (5.38) of (4.20), (4-30), (4-31) is conservative. 

Proof: Let h = BC = 0. Multiplying (5.38) with the integration operator <p T M where <p is smooth, and 
observing that Q sk = ~{Q sk ) T + B^, Q r , = -Qj t + B n ( B B n are defined in (5.15)) leads directly to 

<t> T MJW t - (D\ k <j>) T MF - (IJ rl o) T MG + © I\,)F + <p T (P i © B„ )G = 0. 


The approximation (5.38) is conservative; i.e., it reverses the process of differentiation (second and third 
terms above) and leaves information only at the boundaries (fourth and fifth terms). □ 

5.4.2. Stability. In this Section we will prove the following theorem. 

THEOREM 9. The approximation (5.38) of the problem (4-30), (4-30), (4-31) is both strictly and strongly 
stable if theorem 6 holds and EP r; + /^E < 0. 

Proof: The energy method applied to (5.38) yields 

(5.45) ^(||B'1|3) = GR1 + GR2 + DI + BT E + BT W + BT N + BT S + IT 

where it is assumed that MJ is a norm; the requirements are given in theorem 5. The boundary terms 
BT e + BT\\ BTx + BTs are exactly the same as in the single domain case (see (5.24)), while the 
operator in GR1, GR2. and Dl is replaced by D\ k defined in (5.39). Strict and strong stability of (5.38) 
follows if 


(5.46) IT = W r D m (EP„ + P r ,E)lT < 0. 

Because D > 0, we need EP t/ -b P r; E < 0. □ 

Remark . Because P y] > 0, E < 0 with the first and last r elements in E being constants would satisfy 
condition (5.46). 

6. Numerical experiments. In the calculations below, we have used the fourth- and sixth-order 
schemes reported in [17] in space and a five-stage fourth-order RK scheme [32] in time. 
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Fig. fi.l. Instability due to vanishing wave speed and flux interface conditions. 


6.1. Vanishing wave speed. For problems with a realistic geomet ry, one will frequently encounter zero 
wave speed somewhere in the field due to the variation of the metric coefficients, the variable coefficients, or 
(for nonlinear problems) the solution. This difficulty (see Section 4.3.2) particularly severe in one dimension, 
is exemplified in the calculation of Burger's equation shown in Figure (hi. 

The instability that develops close to zero wave speed when using a penalty on the fluxes at the interfaces 
is evident. With interface conditions applied on the variable instead of the fluxes, the instability disappears. 
Also, if one scales the problem such that V varies between 1 and 3 instead of 0 and 2 one can use flux 
interface conditions without any sign of instabilities. This anomalous behavior associated with a vanishing 
wave speed occurs with other numerical scheme's, and is typically suppressed by adding dissipation (e.g. the 
“Entropy fix'' used with Roe solvers). 

6.2. Error growth due to varying coefficients. Consider the following II) test problem, 

ut + Fr = 0. [.r,.t/]ea />() 

^ ^ u — f, [j\ t/] G ft. / = 0 

on the *20 domain ft = [.r,i/] G [0, 1] x [0. I]. The variables, fluxes and initial data are 


(<>. 2 ) 


Ml 

, F = 

a{j-)n i 

. f = 

I/O 


/>(.(• ) ll-j 
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FlG. 6.2. The error in the growth rate for different transformations. 

The problem (6.1) is 1-periodic in y and has 

■ 0 [£, //] G £1. t > 0 
/. [£, rj\ £ ft, r = 0 

x [0, 1], The problem (6.4) has the same boundary 

conditions as (6.1). 

6.2.1. The energy growth in ID. The energy growth for the ID (y = 0, ?/. r = 0) version of (6.1 )-(6.4) 
with 

(6.5) a=l + ear, b = -i+ex, a = i, $ — y/(l -f e)/(l - f) 

leads to ||«||f — — e||t/|| _ . The growth rate — r/2 corresponds to a single eigenvalue on the real axis in 
the continuous spectrum. Figure 6.2 shows the error in the sixth-order numerical approximation of this 
eigenvalue for different transformations (r = •**(£))• Figure 6.3 shows the convergence (in an Lo sense) of the 
seven eigenvalues with most accurately converged real parts. The convergence rate in both Figures 6.2 and 
6.3 is at least 6. 


(6.3) «i(0,jM) = <*“2(0, ?/,/), 

as boundary conditions in the ^’-direction. 

By introducing a 2D curvilinear mesh we obtain 

(0 4) Ju T + (F)t + (G)n = 

u = 

where F = J£ X F,G = Jr) T F and ft = [£,r/] e [0,1] 
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Seven Most -Accurate Eigenvalues 



FlG. 6.3. The error in the growth rate for varying wave speeds. 

Even though the resolved eigenvalues (and eigenvectors) converge at the theoretical rate (see Figures 6.2- 
6.3), there are unresolved eigenvalues and eigenvectors that can generate difficulties. In Figure 6.4, the least 
resolved eigenvector corresponds to an eigenvalue with a negative real part (-4.6529E-03) significantly more 
to the right of the analytical value (-7.5000E-03) than could he expected by the order of the approximation. 
These unresolved eigenvalues and eigenvectors may generate extra large energy growth, as shown in Figure 
6.5. The growth varies with the initial condition. Note that the extra energy growth for a uniform mesh 
can be present only for varying coefficients because otherwise GR1 — (F, D r V) — (D r F, F) = (). 

6.2.2. The energy growth in 2D. The energy growth for the 21) continuous problems (6.1), (6.4) 
is identically zero with ( = 0 in (6.5); i.e.. the L> norm of the solution remains constant in time. In 
the semi-discrete case, the energy growth is given by (5.45) where GR2 = 1)1 = 0 and the introduction of 
boundary conditions (BTj) and interface conditions (IT) leads to damping. Possible error growth, see (5.1(5), 
is provided by, 

(6.6) UR 1 = -[(r. D ( F) - (D ( C. F)) - [(U. D n G) - (0„I\G)} 

only. For a linear mapping where the metric coefficients are constants (see Figure 6.(5) we obtain I)^F = 
D^\\F — A i Dj: F, D r] (r — D v AoF = A^D V F, which yields GR1 = 0. The error growth is shown in Figure 
6.7. The calculations are fourth-order accurate in time. Note that there is an absolute bound on the error. 

In a nonlinear mapping (see Figure 6.8) the truncation errors in the metric calculation, and consequently 
also in the calculation of the fluxes, leads to GRl^ 0, which in turn can generate error growth (see Figure 



SYSTEM STABILITY: EIGENVECTORS 


Physical, First, and Ra(max) 



FlO. 6.4. Eigenvectors related to two most and least resolved eigenvalues. 


6.9). These calculations also are fourth-order accurate in time. Note the enormous time scale in Figures 6.7, 
and 6.9. 


6.3. Navier-Stokes calculations. We consider here a ID viscous shock propagating in accordance 
with a Mach number of 2.0 and a Reynolds number 150 over a 2D domain. The exact solution of the 
Navier-Stokes equation for this case can be found in [33]. At the artificial boundaries, including the circular 
region in the middle, we impose flux boundary conditions by using the penalty formulation on the fluxes 
with exact data from the analytical solution. At the interfaces we impose interface conditions by using the 
penalty formulation on the variables. 

In Figure 6.10, the density and grid for the propagating shock is shown. The shock travels from the lower 
left corner to the upper right corner and has almost passed out of the computational domain that consists of 
12 blocks. The sixth order scheme and 24 gridpoints were used in each sub-domain. The local density errors 
are shown in Figure 6.11. The grid refinement study in Table 6.3 indicate between fifth- and sixth-order 
accuracy in an Lo norm, consistent with the theory in [34], [35], since we have fifth order accuracy at the 
boundaries and interfaces (see (3.4)) and relatively coarse grids. 

7. Summary and conclusions. We have analyzed boundary and interface conditions for high or- 
der finite difference methods applied to multidimensional linear problems in curvilinear coordinates. The 
investigation focused on the effect of variable coefficients. 

A problem with a norm as a function of the Jacobian was analyzed. Boundary and interface conditions 
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FiG. (>.5. Extra growth due to unresolved features and initial conditions. 


Table 6.1 

Twelve subdomains , sixth- order explicit; CFL = 0.3. 


Wave speed 

49/65 

65/97 

97/129 

129/19:5 

-0.25 

-4.610 

-4.640 

-4.722 

-4.722 

0.00 

-5.115 

-4.986 

-4.538 

-4.657 

0.25 

-5,155 

-5.253 

-5.179 

-4.952 

0.50 

-5.331 

-5.401 

-5.467 

-5.327 

0.75 

-5.523 

-5.514 

-5.590 

-5.565 

1 .00 

-5.635 

-5.622 

-5.659 

-5.719 

average 

-5.22* 

-5.236 

-5.193 

-5.196 


in both flux and variable formulations have been investigated. Flux boundary conditions lead to energy 
estimates whereas flux interface conditions lead to difficulties if the wave speed approaches zero. 

A new and simplified so called Q-formulation of the penalty method was derived at interfaces. The 
Q-formulation simplifies and extends the formulation and implementation of derivative conditions in both 
one and two dimensions at interfaces. 

It was shown that varying coefficients can cause unbounded error growth via the truncation errors 
even though the boundary and interface conditions are implemented in a stable and dissipative way. The 
error growth may be large due to unresolved features in the solution. Numerical calculations confirmed the 





Fig. 6.6. A 4 block mesh, linear mapping. 


theoretical conclusions. 
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FlG. 6.10. Propagating viscous shock. 
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Fig. fi.ll. Local error levels for propagating viscous shock. 
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